Bohmian-trajectory analysis of high-order harmonic generation: 
Probing the local dynamics 
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In this article, we perform detailed studies of high-order harmonic generation (HHG) employing 
Bohmian trajectories, with particular emphasis on which spatial regions are the most relevant for this 
phenomenon. These trajectories, directly extracted from the time-dependent Schrodinger equation, 
play the role of "tracer particles" , thus being very helpful for mapping the flow of the time-dependent 
wavefunction. In particular, they allow us to revisit several issues, such as the use of the dipole 
length or the dipole acceleration in the computation of harmonic spectra, the influence of numerical 
filters on the fiow of the wavefunction in real time, or the role of the innermost and outermost core 
regions and of ionization in the spectra. Our results show that the Bohmian trajectory starting 
at the innermost part of the core is the most relevant in order to reproduce typical HHG spectra, 
regardless of the pulse shape and of the potential range employed. This provides additional insight on 
the success of strong-fleld models that consider the atom to be a source term, such as the strong-fleld 
approximation. A quantitative agreement between the spectra obtained from Bohmian-trajectory 
computations and the time-dependent Schrodinger equation is also obtained. Furthermore, we also 
establish a parallel between the use of the dipole acceleration and the Banning flher, by analyzing 
how, in both cases, Bohmian trajectories away from the core are removed. 

PACS numbers: 32.80.Rm, 42.50.Hz 



I. INTRODUCTION 

High-order harmonic generation (HHG) is a weU- 
known phenomenon that occurs when matter interacts 
with a strong laser field, of intensities of the order of 
or higher than / ~ 10^^ W/cm^. Since the late 1980s, 
this phenomenon has been investigated both theoretically 
and experimentally. It has attracted a great deal of at- 
tention due to the myriad of possible applications, such 
as table-top, extreme ultraviolet (XUV) sources [IJ, high- 
frequency light pulses of attosecond duration (for a 
review, see Ref. Q), and, in recent years, it has also been 
identified as a powerful tool for the attosecond imaging 
of dynamic processes in matter 5-9]. These applications 
have been made possible due to a very intuitive phys- 
ical interpretation of HHG in terms of electron orbits, 
in which an electron, under the influence of the external 
laser field, reaches the continuum by tunnel ionization, 
is accelerated by the field, and, at a subsequent time, 
recollides with its parent ion TO]. Upon recollision, the 
electron may recombine with a bound state, thus emit- 
ting high-frequency radiation. This physical picture has 
been put across in the context of the three-step model 
(TSM). It has been explored both in classical models 
[lO| and in quantum-mechanical a ppr oaches, such as the 
strong- field approximation (SFA) Furthermore, the 
recollision picture has also been inferred from purely nu- 
merical models, such as the numerical solution of the 
time-dependent Schrodinger equation (TDSE), by time- 
frequency analysis (see, e.g., |12Hla |). 

The above-stated picture is consistent with the fact 
that, in the modeling of HHG, to first approximation, it 



suffices to assume that the core is static and to treat it as 
a source term, located at the origin r = of the coordi- 
nate system. This is the key assumption behind the SFA, 
which approximates the continuum by field-dressed plane 
waves. On the one hand, this simplifies the problem con- 
siderably and allows for intuitive, semi- analytical models 
to be developed. On the other hand, it does not allow 
for quantitative predictions of harmonic spectra and un- 
derplays the role of the core. In fact, even for a simple, 
one-electron system, such as the hydrogen atom, pulling 
out one electron will alter the charge distribution in the 
core considerably. On a more general note, it has been 
recently shown that core polarization and excitation Q , 
as well as the infiuence of the Coulomb potential, are ex- 
pected toplay an important role in strong-field phenom- 
ena [H, [l3| . Specific examples are the so-called low en- 
ergy structure [Tsl and interference patterns [l^ in ATI. 

Ab initio methods, such as the TDSE, though, provide 
the full time-dependent wavefunction, but do not give a 
clear space-time picture in terms of electron orbits. Typi- 
cally, this picture is extracted by using windowed Fourier 
transforms and constructing time-frequency maps. Oth- 
erwise, the main focus is on how to minimize the numer- 
ical background and spurious refiections in the spectra. 
From this perspective, it is well known that using the 
expectation value of the dipole acceleration instead of 
that of the dipole length leads to higher- quality spec- 
tra [13, Physically, it has been argued that, since 
the light emitted by a moving charge, i.e., the electron, 
is proportional to its dipole acceleration, it would make 
more sense to use this quantity as the observable from 
which HHG spectra are calculated. On a more techni- 
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cal level, the dipole acceleration probes regions near the 
core, while the dipole length emphasizes large distances. 
This implies that the latter observable is far more sen- 
sitive with regard to reflections at the edges of the box 
employed in the numerical integration. This can be re- 
duced using absorbing boundaries and a Hanning filter, 
but obviously one also expects that these devices will al- 
ter the system dynamics. 

Apart from the aforementioned approaches, there exist 
other theoretical methods, which allow a space-time pic- 
ture of HHG in terms of electron orbits and still involve a 
lesser degree of approximation than the SFA. One of such 
methods is Bohmian mechanics [12, ■ In strong- field 
physics, so far, Bohmian mechanics has only been em- 
ployed very recently in a handful of publications [23 - l30| . 
It is, however, a powerful tool in order to assess the in- 
terplay between the residual potential and the laser field. 
Within an analytic working scheme [3l| . Bohmian trajec- 
tories are directly obtained from the full solution of the 
TDSE. Hence they do not come from any sort of phys- 
ical approximation, but they contain all the quantum 
phase information necessary for the full description of 
the problem. Because of this reason, Bohmian trajecto- 
ries are employed in many areas of physics and chemistry 
as "tracer particles" [32, i33| , which provide an insightful 
hydrodynamic-like picture of the probability-density flow 
associated with the quantum mechanical wavefunction 
(indeed, this picture is pretty much related to Madelung's 
former hydrodynamical reformulation of quantum me- 
chanics |34i]). 

Recently, we have computed HHG spectra employing 
Bohmian trajectories [30[, and compared such results 
with those obtained from the numerical solution of the 
one-dimensional TDSE. We found that the trajectory lo- 
cated at the innermost part of the core [i.e., starting at 
x(0) — 0] best reproduces the HHG spectra obtained 
from the dipole acceleration. Both spectra were quali- 
tatively very similar, apart from an overall difference in 
intensity. This provided additional insight on the role 
played by specific spatial regions in HHG. First, since 
the innermost region contributes the most, it is legitimate 
to assume that, to first approximation, the core may be 
represented by a source term at the origin, a: = 0, as con- 
sidered in the SFA. Second, using time-frequency maps, 
we have shown that a spatially confined Bohmian trajec- 
tory may contain both bound and continuum dynamics. 
In fact, the specific time-frequency profile obtained from 
the above-mentioned Bohmian trajectory was associated 
with that of an ensemble of classical trajectories under- 
going an unbound motion. In contrast, the motion of 
bound classical trajectories may be associated with spe- 
cific, confined regions in space. A still open question, 
though, is what specific regions in space are necessary to 
obtain a quantitative agreement between the Bohmian- 
trajectory and the full TDSE computations. 

In the present work we discuss a number of issues that 
will give a robust idea about why Bohmian mechanics 
is an excellent tool to probe the local dynamics in the 



HHG phenomenon. To that end, we focus the discussion 
on the dipole length vs dipole acceleration issue, and pro- 
vide an alternative insight on the interplay between the 
electron and the core dynamics in HHG. Accordingly, we 
have organized this work as follows. In Sec. |lT] we pro- 
vide the necessary theoretical background to understand 
our subsequent results, which includes the model con- 
sidered, some remarks about Bohmian mechanics, the 
choice of initial conditions, and the propagation meth- 
ods used here. The results obtained are presented in 
Sec, mil More specifically, in Sec. IIIIBj we commence by 
making a qualitative assessment of the different contri- 
butions to the spectra for specific trajectories, near the 
core and near the boundary. Thereafter, in Sec. IIII CI 
we present ensemble computations, in which we identify 
which regions of the core lead to a quantitative agree- 
ment between the Bohmian-trajectory computation and 
the TDSE, for both length and acceleration. Finally, in 
Sec. IIII D1 we assess the influence of Hanning filters on the 
electron dynamics, far away and in the core region. Here 
we show that Bohmian mechanics allows a reconstruction 
of the equivalent dynamics to having a filter, something 
which is not possible at all with the wavefunction. Note 
that once the expectation value (either the length or the 
acceleration) is computed to obtain the corresponding 
spectrum, all the information about the local configura- 
tion space dynamics is irreversibly lost (in other words, 
one cannot obtain the wavefunction or the probability 
density from such an expectation value). Closing this 
work, the main conclusions and results are summarized 
in Sec. |lVl 

II. THEORY 

A. The model 

For the sake of clarity, we have developed our analysis 
in one dimension, which suffices to illustrate the time- 
dependent electron dynamics, and yet contains the es- 
sential physical elements for linear polarization. Making 
use of the dipole approximation and the length gauge, 
the electron dynamics is then described by the time- 
dependent Hamiltonian (in atomic units) 

H = -\^^^^V{x)-xE^F{t), (1) 

where EoF{t) denotes the driving field and V{x) the 
binding potential (note that the minus sign in the last 
term arises from the electron charge: in atomic units, 
e = -l). 

The quantities that we are interested in here are the 
dipole length, 

x{t)^{^\x\^), (2) 

which, in atomic units, is equal to the dipole moment, 
but with opposite sign, since /i(t) = ex{t), and the dipole 



3 



acceleration, 



-a{t) = -{^\V'ixm, 



(3) 



where V'{x) = dV{x)/dx is defined as the acceleration 
operator. In these expressions, "^f^x^t) denotes the wave- 
function solution of the TDSE, 



ih —^^H^{x,t), 



(4) 



with H being the Hamiltonian operator ([T]). Following 
the standard procedure, from these quantities one ob- 
tains the corresponding power spectra, 



lie) 



9{t)e- 



(5) 



where g{t) generically denotes either x{t) or a{t) [or a 
Bohmian trajectory, x{t), as seen below]. 



B. Bohmian trajectories 

In order to determine the Bohmian trajectories, first 
the above wavefunction is recast in polar form. 



(6) 



where p is the probability density and S is the real- valued 
phase. After substituting this form into the TDSE ^ 
and rearranging terms, one obtains a set of two real- 
valued differential equations, which account for the evo- 
lution of these two quantities. One is the continuity equa- 
tion. 



dp 
dt 



V - J = 0, 



(7) 



where J = pS/S/m is the usual quantum probability cur- 
rent density. The other equation is a quantum Hamilton- 
Jacobi equation. 



dS , (V5)^ 



2m 



V + Q = 0, 



(8) 



where 



Q{x,t) = 



'2m pi/2 



4m 



P 



(9) 

is the so-called quantum potential and S is the equivalent 
of the classical action 

As happens with the classical analog of Eq. ([8]), one 
can postulate [l^, the existence of trajectories (char- 
acteristics) that are solutions to this equation. As in the 
classical case, these trajectories are perpendicular to sur- 
faces of constant phase and therefore the corresponding 
equation of motion is proportional to VS. Alternatively, 
the same can also be straightforwardly stated if Eq. ([7]) 



is regarded as describing a (quantum) transport problem 
[121, with a local velocity field being defined as v — J / p. 
In either case one obtains the equation of motion 



J 

P 



h 

2im 



(10) 



which, after integration, renders the corresponding 
Bohmian trajectories, i.e., the paths along which the 
probability flows. 

Once the Bohmian trajectories are synthesized, two 
types of calculations have been carried out with them. 
The first type consists in the direct evaluation of their 
individual spectrum in order to determine how each tra- 
jectory contributes to the expectation values and ([3]). 
That is, we have studied the power spectrum ([S]) for a 
series of individual trajectories, substituting g{t) by the 
time-dependent trajectory x(t) obtained after integrating 
Eq. (|10l) with the corresponding initial condition, x(0). 
In a second step, we have stressed the ensemble feature 
of Bohmian mechanics, analyzing the emergence of the 
usual dipole length and acceleration as the number of 
Bohmian trajectories was increased [and bearing in mind 
that their initial positions distributed randomly accord- 
ing to ^(x, 0)]. In this case, we have taken into account 
the fact that the Bohmian version of the expectation val- 
ues ^ and dS]) are given by 



1 ^ 



1 ^ dV{x) 



N ^ dx 

i=0 



(11) 

(12) 



=Xi{t) 



respectively, which follow a usual statistical approach. In 
these expressions, Xi{t) denotes the ith Bohmian trajec- 
tory from an ensemble and N is the total number of these 
trajectories in each numerical experiment performed. 



C. Initial conditions and dynamics 

As the initial condition for the wavefunction, we have 
considered the ground state of the bare or field-free 
Hamiltonian, 



V{x) 



(13) 



This eigenstate has been numerically obtained by means 
of the imaginary time propagation method 35[. Then, 
the exact time-propagation of the wavefunction accord- 
ing to the TDSE (j4]) has been carried out by combining 
the split-operator technique ^^'] with the fast Fourier 
transform (EFT) technique [37|. More specifically, the 
latter has been used to achieve a high performance in 
the evaluation of the action of the kinetic operator on 
the wavefunction. We have also taken advantage of this 
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decomposition of the wavefunction to recast purposely 
the wavefunction as 



ipn x/h 



(14) 



Here, the c„(t) coefficients denote the Fourier weights 
(in general, complex valued) associated with each one of 
the plane-wave components with momentum pn of the 
decomposition at time t. 

To ensure that all the relevant dynamics are incorpo- 
rated for the parameter range of interest, we have set 
the box boundaries located far enough from the core re- 
gion (at |a;| — 150 a.u.). Furthermore, special care has 
been taken in order to avoid reflections and spurious ef- 
fects near the box edges. First, we have employed a 
mask function in the form of a cos^^^-function, which 
becomes active near the boundaries. This function is 
smooth enough, but still capable to absorbing with a high 
efficiency, thus avoiding nonphysical reflections. We have 
tested this fact by considering a wide range of box sizes. 
Second, the maximum/ minimum value of the total po- 
tential function has been truncated in order to avoid also 
nonphysical accelerations towards the box edges. In this 
regard, the size of the box was chosen such that this trun- 
cation takes place close to the region where the absorber 
becomes active. 

Additionally, to avoid effects associated with the loss 
of probability due to the presence of the absorber, which 
becomes more prominent for higher values of the exter- 
nal field intensity, in the calculations special care was 
taken in properly renormalizing the wavefunction. Fur- 
thermore, unless otherwise stated (see Sec. IIII D|) . no fil- 
ter functions (e.g., Hanning windows) were used in order 
to avoid influencing the long-range dynamics and there- 
fore the topology of the Bohmian trajectories. 

Regarding the Bohmian equation of motion (fTO|) . it has 
been integrated "on the fly" , substituting the value of the 
wavefunction 4'(a;,t) at each time into the last term of 
this equation. Now, given the general form (fT4)) . Eq. ^TU\\ 
can be recast as 

1 J2mnP"\^"ii't)\\^ri{t)\cOs[ApnmX/h+ Anmit)] 
X — ' 

I]mnlcm(OI|cn(i)|cos[Ap„„ia;/?i + A„„j(i)] 

where the Cn(t) coefficients have been expressed in po- 
lar form, i.e., c„(i) — |c„(t)|e*''"(*\ with |c„(<)| and 5„(t) 
denoting its (time-dependent) modulus and phase, re- 
spectively. Moreover, we have defined Apnm = Pn — Pm 
and Anmit) = <5„(t) — Sm{t). In this way, since the c„(i) 
are known at each time through the FFT operation, the 
Bohmian trajectories x{t) are obtained straightforwardly 
by integrating the analytical form (jlSp , which is an ordi- 
nary differential equation. In particular, in this work the 
integration has been performed by means of the Euler 
method, which has been proven to be accurate enough. 

Unless otherwise stated, the initial conditions a;(0) 
considered for the Bohmian trajectories have been ran- 
domly chosen to mimic the initial probability distribu- 
tion, p{x,0) = |\['(a;,0)p, within the interval [—Xc,Xc], 



with Xc = 4.102 a.u. This ensures that the time-evolution 
of most of the probability density will be well monitored, 
as the integrated probability for |x| > Xc is only 0.151% 
of the total probability. 



III. RESULTS AND DISCUSSION 

A. Driving fields and binding potentials 

In the results that follow, we will mainly assume that 
the atomic potential, V{x), is given by the soft-core in- 
teraction 



1 



(16) 



where a is the impact parameter. Far away from the 
core, this potential decays as x^^, behaving similarly to 
a Coulomb-type, long-range interaction. Equation (jl6p 
gives dV{x)/dx = x/{x^ + a^)^/^ for the dipole accelera- 
tion operator in Eq. (jS]), and the acceleration along 
the Bohmian trajectories reads as 



N 



iV^[l+xf(t)]3/2- 



(17) 



Here we have chosen the impact parameter to be a = 1, 
which yields a ground-state energy eo = —0.66995 a.u. 
for the field-free Hamiltonian (IT^ . 

In Sec. IIIIB 21 the long-range of the potential V^aT(^) 
is eliminated by re-writing it as 



VUx) = - 



1 



■fix), 



\/x^ + a 

where the truncating function is described as 



(18) 



1, 



f{x) = { cos"^ 
0, 



7 ( TT l^l - ap 
2 L — ag 



\x\ < ao 
ao<\x\<L . (19) 
\x\ > L 



The explicit values of ag and L are chosen in two different 
ways, which will be later on specified in that Section (see 
Table U for details) . A similar truncation was formerly 
used in Ref. [S^ to assess the influence of the potential 
tail on resonant HHG enhancements. 

Regarding the external laser field, F{t) — f{t) sm{ujot), 
two types of envelopes, f(t), have been considered, 
namely a flat-top trapezoidal pulse, 

f (t/ro), 0<i<To„ 

fit) ^ I I, Ton < i < Toff , (20) 

[ 1 - (i - Toff)/ro„, Toff < < < Tf 

and a sin^-function, 

/,2(i) =sin2(7ri/T/). (21) 
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In both cases, wq = 0.057 a.u. and tq = 2tt/ujo denote 
the frequency and period of the field, respectively. More 
specifically, in ([20l) the time required for the external field 
to reach its maximum amplitude is Ton — 2.25 tq, while 
Toff = 2.25 To + Ntq stands for the duration of the field 
with maximum amplitude, and tj = Tog + 2.25 tq refers 
to the time at which the field is switched off (and our 
quantum simulations are completed); to ensure a high 
monochromaticity in the flat-top pulse, we have chosen 
TV = 10 (t/ = 14.5to). The propagation described in 
Sec. Ill CI has been taken up to t — Tf. 

B. Single-trajectory analysis 

1. Soft-core potential and flat-top pulse 

We will first revisit the issue of the dipole length vs 
acceleration for the soft-core potential ([T6l) employing 
Bohmian dynamics. The dipole length over-emphasizes 
regions near the box edges due to its linearity with the 
distance. Hence, even if the probability density far from 
the core is negligible, x can reach relatively high values, 
thus compensating this fact and enhancing the associ- 
ated features in the corresponding spectra. In contrast, 
the dipole acceleration focuses on the core region. Specifi- 
cally for the potential ([T5| . it increases up to a: = ±l/-\/2, 
and then decreases to zero again very rapidly, like x~^. 
Hence, any contribution from portions of the probability 
density outside this region is strongly suppressed f2^. 

In Fig. [U we show the time-dependent dipole length 
[panel (a)] and acceleration [panel (b)] obtained from 
Eqs. ^ and ([3]), respectively, and computed using the 
TDSE. Moreover, the corresponding harmonic spectra 
are provided in the lower panels. As seen in Fig. [Ija), 
as the field strength Eq increases, the amplitude of the 
oscillations for the dipole length becomes comparatively 
large. Furthermore, at the end of the pulse, the expecta- 
tion value of the dipole operator is non- vanishing for the 
larger driving-field intensity. This can be traced to irre- 
versible ionization and the contributions far from the core 
region, which are over-estimated in this case. The dipole 
acceleration [see Fig.[Ub)], however, behaves in a rather 
different way. First, there is a much less pronounced 
increase in the amplitude of a{t) with the driving-field 
strength. Comparing the two graphs, we note that their 
main difference lies in the substructure observed for a(t). 
Second, the acceleration vanishes at the end of the pulse. 
These differences stem from the fact that a{t) only de- 
pends on the part of the wavefunction where the electron 
undergoes the strongest interaction with the core, as dis- 
cussed above. 

When comparing the spectra for the dipole length 
and the acceleration for such intensities, however [see 
Figs. [T{c) and (d)], one readily notices that the spec- 
tra obtained from the dipole acceleration exhibit a large 
plateau consisting of high-order harmonics extending up 
to a cutoff energy of Wc ~ |eo| -I- 3.17Up, where Up = 



0.12 




10 20 30 40 50 60 70 80 10 20 30 40 50 60 70 80 



Harmonic order (oj/coq) Harmonic order (co/wo) 

FIG. 1: (Color online) Dipole length (a) and dipole acceler- 
ation (b) computed from the TDSE for a model consisting 
of a soft-core atomic potential and a flat-top laser pulse with 
frequency ujq = 0.057 a.u., and intensities Eq = 0.05 a.u. 
(dashed red line) and Eq — 0.075 a.u. (black solid line). The 
corresponding power spectra are displayed in panels (c) and 
(d), respectively. The cutoff frequencies in the cases consid- 
ered are uj — 22.5cjo for Eq = 0.05 a.u. and lj = 35.8a;o for 
Eo = 0.75 a.u.. 



Eq/4ujq is the ponderomotive energy. These features, 
shown in Fig. [IJd), are in full agreement with the pre- 
dictions of the TSM. In the harmonic spectra computed 
from the dipole length, however, these features are not so 
easily identifiable. In particular, for the higher driving- 
field intensity, the harmonic peaks are obfuscated by a 
large background, with no visible plateau. Overall, the 
HHG spectra in Fig. [Ijc) are monotonically decreasing. 
Nonetheless, for moderate intensities of the external field 
(see curve for Eq — 0.05 a.u.) the cutoff could be as- 
signed to the dramatic decrease in the intensity observed 
around uj « 30wo- 

We will now employ Bohmian trajectories in order 
to analyze in detail why the spectra obtained from the 
dipole acceleration and from the dipole length are so dif- 
ferent (in particular, why the cutoff frequency does not 
appear clearly in the spectra obtained from the dipole 
length). In Figs. Ufa) and (b) two sets of Bohmian tra- 
jectories are displayed for the same driving-field intensi- 
ties as in Fig. [TJ As Eq increases some of the trajectories 
start to oscillate with increasingly larger amplitudes and 
eventually leave the core region [see panel (b)]. Actually, 
it is this behavior with Eq what leads to the eventual de- 
pendence on the boundaries of the problem. That is, as 
Eq increases, the trajectories undergo larger excursions 
out of the core and therefore, when their average (fTTj) is 
computed, one will obtain results highly dependent on 
the boundary conditions chosen to solve the TDSE. Note 
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10 20 30 40 50 60 70 80 10 20 30 40 50 SO 70 80 



Harmonic order Ito/tan) Harmonic order (to/o)n) 

FIG. 2: (Color online) Top: Bohmian trajectories obtained 
with a flat-top pulse with intensity: (a) Eo — 0.05 a.u. and 
(b) Eo — 0.075 a.u. Bottom: Power spectra computed from 
individual trajectories Xi{t) selected from among the sets at 
the top panels, (c) Trajectories starting at a::(0) = a.u. (black 
solid line) and 2;(0) — —3 a.u. (red dashed line), for _Eo = 
0.05 a.u. (d) Trajectories starting at a::(0) — a.u. (black 
solid line), 2;(0) = 1.8 a.u. (blue dotted line) and a::(0) = 
—3 a.u. (red dashed line), for Eo = 0.075 a.u. Colors/line- 
styles are in correspondence with the trajectories displayed 
at the respective top panels. 



that this average formally corresponds to the expectation 
value of the dipole moment. On the contrary, if one only 
focusses on those trajectories that always keep oscillating 
around the core, the average behavior will resemble that 
of the dipole acceleration. 

The interpretation provided above is supported by the 
harmonic spectra obtained from single Bohmian trajecto- 
ries Xi{t), depicted in the lower panels of the figure. As 
an overall feature, a large plateau followed by a sharp 
cutoff is always obtained if the trajectory starting at 
x(0) = is taken, regardless of the driving-field inten- 
sity. This trajectory corresponds to the innermost part 
of the wavefunction, which is strongly localized at the 
core. This spatial region resembles that spanned by the 
dipole acceleration. In contrast, the spectra from tra- 
jectories whose initial positions lie at a few atomic units 
of distance, i.e., at the outermost regions of the wave- 
function, will be strongly infiuenced by the driving-field 
intensity. For moderate driving-field intensity [Fig. [2jc)] 
the spectra from both innermost and outer trajectories 
are quite similar in shape, apart from an overall intensity 
difference. This is caused by the fact that both trajec- 
tories perform an oscillatory motion near the core. If, 
however, the intensity increases, the spectra of the inner- 
most and of the outermost trajectories differ markedly, as 
exemplified in Fig.[21[d). In the figure, one observes that. 



if the trajectories are launched from a;(0) = —3 a.u., or 
even a;(0) = 1.8 a.u., the corresponding spectra will only 
consist of the fundamental and of a uniform background. 
For these sets of initial conditions, the confined oscillat- 
ing motion (around the corresponding initial position) is 
lost. 

Apart from that, the further from a;(0) =0 the initial 
position of a trajectory is, the higher the overall inten- 
sity in the spectra will be. Thus, even if the number of 
trajectories that leave the core region is relatively small, 
the fact that their spectra are several orders of magni- 
tude more intense than those obtained from the inner- 
most trajectories gives rise to a masking of the latter 
contributions in the dipole-length spectrum. These find- 
ings are in agreement with the observations in Ref. '2^, 
where a background in the HHG spectra was attributed 
to irreversible ionization at the end of the pulse, and in 
Ref. [2l|, where this background was attributed to the 
probability density near the edges of the integration box. 



2. Short-range potentials and flat-top pulse 

In order to investigate the influence of confinement fur- 
ther and whether the above-discussed picture remains 
valid across all model potentials, we will consider the 
truncated short-range potentials given by Eq. ([T5)) "385. 
In this section, the parameters oq and L in Eq. (jl8|) are 
chosen in two distinct ways. For the truncated poten- 
tial V^j^\x), the long-range tail is eliminated, but the 
core region is kept practically unaltered, i.e., the field- 
free eigenenergies obtained for the untruncated potential 
V^^ and V^^\x) are very close, and many weakly bound 
states are present. In contrast, the truncated potential 
(x) is so short ranged that it can only support only 
two strongly bound states. Throughout, we keep the 
ground-state energy roughly the same. A list with the 
eigenvalue spectrum for these three potentials is given in 
Table H 

These potentials, together with the soft-core model 
(ITBl) . are illustrated in Fig. [31 Panels (a) and (b) show 
the field-free potential V{x) and the effective potential 
barrier Vcti{x) = V{x) — xEq when the driving laser field 
reaches its peak value Eq — 0.075 a.u.. This higher in- 
tensity will be employed throughout this section, as it 

causes more irreversible ionization. For Vf^\x), the po- 
tential well (core region) is essentially the same, while the 
long-range tail is significantly decreased. In the presence 
of the field, however, Vcff{x) is practically not affected 
by the truncation. For the second truncated potential 
Vf). (x), the width of the potential well has decreased re- 
markably as well as the potential range. This leads to a 
much higher ionization barrier Vcsix). 

In Fig. |4j we display the dipole length and the acceler- 
ation computed with the TDSE for v}^\x) and v}^\x), 
together with the corresponding power spectra [upper 
and lower panels, respectively]. In Fig. Hl^a) we observe 
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TABLE I: Eigenvalues for the untruncated soft-core (|16|l and 
two truncated soft-core (|18l) potentials. We shall refer to these 



models as first truncated, or v}^\x), for which ao = 5.0 and 
L — 50, second truncated, or V}^\x), for which ao = 1.0 and 
L = 7.55. Note that, in principle, the number of eigenstates 
supported by the untruncated potential approaches infinity; 
in our calculations, though, we obtain a finite number of them 
because of the boundaries of the grid we are using to solve the 
TDSE. All quantities are given in a.u. 



n 


untruncated 


truncated 1 


truncated 2 





-0.66995 


-0.66995 


-0.65651 


1 


-0.27508 


-0.27503 


-0.17864 


2 


-0.15158 


-0.15059 




3 


-0.09276 


-0.08714 




4 


-0.06358 


-0.05013 




5 


-0.04552 


-0.02390 




6 


-0.03462 


-0.00754 




14 


-0.00826 






15 


-0.00707 






16 


-0.00670 







that x{t) undergoes longer excursions with the first trun- 
cated potential V^^\x) than even with the untruncated 
case. Furthermore, there is a marked increase in the 
dipole-length amplitude with regard to the time, which 
is absent in Fig. [Tfa) . Both effects suggest an increase 
in ionization for vj^^^ {x). This behavior can be explained 
as follows. First, both the untruncated potential and 
v//^ (x) support many loosely bound states. In the pres- 
ence of the field, these states will be strongly coupled 
to the continuum. Furthermore, for V^^\x) there will 
be no attractive long-range tail pulling the ionized parts 
of the electronic wavepacket towards the core region. 
This enhanced ionization leads to a large background in 
the corresponding harmonic spectra [see dashed line in 

Fig. m^c)]. For the truncated potential V^^ {x), in con- 
trast, x{t) essentially oscillates within the core region. 
This happens because this potential only contains two 
bound states, the excited one being also quite far from 
the continuum. Hence, the ionization barrier remains rel- 
atively high and, in the presence of the field, there is less 
coupling with the continuum. This outweighs the absence 
of a long-range tail mentioned above. This suppression 
of ionization considerably reduces the background in the 
harmonic spectra and the overall harmonic intensity [see 
solid line in Fig.|4l[c)]. 

In this regard, one should notice that the expectation 
values of dipole acceleration obtained for both truncated 
potentials, displayed in Fig. UJ^b), are very similar. This 
is expected, as the acceleration probes the spatial regions 
close to the core and the main differences are due either 
to the high-lying excited states, whose spatial extension 
is large, or to the different potential barriers. For the very 
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FIG. 3: (Color online) (a) Comparison between the soft- 
core model (|16|l (black solid line) and the truncated soft-core 
model (I18p for two different values of the parameters ao and 
L. The blue dashed line denotes the truncated potential with 
ao = 5.0 and L = 50, while the red dotted line stands for 
ao = 1.0 and L = 7.55. In panel (a) and (b) we display the 
same potentials, but acted by the laser field at the maximum 
of its amplitude, which is taken as Eo — 0.075 a.u.. 
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FIG. 4: (Color online) Dipole length (a) and dipole accel- 
eration (b) for the truncated soft-core model (|18|l with two 
sets of parameters: ao = 5.0 and L = 50 (red dashed line), 
and ao — 1.0 and L = 7.55 (black solid line), computed from 
the TDSE. The pulse considered was a flat-top one, as given 
above, with intensity Eo — 0.075 a.u. and the same frequency 
as in the previous frequencies. The corresponding power spec- 
tra are displayed in panels (c) and (d), respectively. For the 
potentials V}^\x) and V^^^x), the predicted cutoff energy 
lies at cj = 35.8aJo and u — 35.6a;o, respectively. 



same reason, the harmonic spectra exhibit a much clearer 
cutoff, near u — 35u!q and similar shapes [see Fig. Sfd)]. 
Furthermore, the discrepancy in the overall harmonic in- 
tensities is not as extreme as in Fig. HJc). Still, due to 
the suppression of ionization, the spectrum obtained with 
{x) is a few orders of magnitude weaker. 
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More detailed information about the ionization dy- 
namics is provided by the outcome of the Bohmian- 
trajectory computation, exhibited in Figs. EJa) and (b). 

For the short-range potential v}^^ (x) , the largest amount 
of ionization is essentially produced in the central part 
of the pulse, as seen in panel (a) , while in the long-range 
soft-core model (|16p we found that ionization was pro- 
duced in a more homogeneous way along the pulse du- 
ration. This is related to v}^\x) being short ranged: 
once the trajectories are outside the core region, there is 
no potential tail attracting them towards it. This leads 
to an increase in the outwards probability flow. On the 
other hand, for the second short-range potential Vj^ {x), 
suppression of ionization due to a higher barrier and to 
the lack of loosely bound excited states plays a more im- 
portant role. Consequently, most trajectories oscillate 
near the core region. Nevertheless, given the intensity of 
the laser held, eventually some trajectories escape lead- 
ing to irreversible ionization. For the same reason, i.e., 
the shorter range of the potential, they will remain in the 
continuum in spite of the laser field amplitude becoming 
very small near the pulse turn off, as shown in panel (b) . 
Similarly to what has been observed for the long-range 
potential, the central trajectory, starting at x(0) = 0, is 
the one that best reproduces the behavior displayed by 
the dipole acceleration, while marginal trajectories ren- 
der spectra which decrease monotonically, except for a 
single peak at the field frequency, uj = ujq [see panels (c) 
and (d)]. 
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FIG. 5: (Color online) Top: Bohmian trajectories for the 
truncated soft-core model (jlSp with two sets of parameters: 
(a) ao — 5.0 and L — 50, and (b) ao = 1.0 and L = 7.55, in 
a flat-top pulse with intensity iJo = 0.075 a.u. and the same 
frequency as in the previous flgures. Bottom: Power spectra 
corresponding to some selected trajectories from among the 
sets, (c) Trajectories starting at a;(0) = a.u. (black solid 
line) and a;(0) = —3 a.u. (red dashed line), (d) Trajecto- 
ries starting at a;(0) = a.u. (black solid line), a;(0) = -1-3 
a.u. (dotted blue line) and a;(0) = —3 a.u. (red dashed line). 
Colors/line-styles are in correspondence with the trajectories 
displayed at the respective top panels. 



3. Soft-core potential and sinusoidal pulse 

Another interesting question that arises here is that of 
the role played by the overall shape of the laser pulse. 
Up to now we have considered a flat-top pulse, which to 
some extent may be viewed as a reasonable approxima- 
tion to a monochromatic wave. One could argue, how- 
ever, that the pulses employed in experiments are far 
from monochromatic, and the pulse shape influences the 
dynamics. For that reason, we will perform a Bohmian 
trajectory analysis for the long-range potential model 
(1161) . but assuming now that the driving field is a pulse 
whose envelope is given by Eq. (Hi]). 

The outcome of our computations is shown in Fig. |51 
for the two peak strengths employed in the remaining 
sections of this work. In the left panels we display the 
dipole acceleration and the corresponding power spec- 
tra obtained from the TDSE, which will be employed 
as benchmarks. In general, the HHG spectra follow the 
cutoff law |eo| + 3.17t/p. However, the resolution of the 
harmonic peaks is much worse than that obtained for the 
fiat-top pulse. Indeed, only in the cutoff region it is pos- 
sible too identify the harmonic peaks. This is expected as 
sharp harmonic peaks stem from the field periodicity: a 
trapezoidal pulse contains many identical cycles, whereas 
a sinusoidal pulse does not. The remaining panels of the 



figure exhibit the Bohmian trajectories obtained for this 
kind of pulse [panels (b) and (c)], with the correspond- 
ing spectra [panels (e) and (f)]. For sinusoidal driving 
pulses, we find the same pattern described above, namely 
that there is a good qualitative agreement between the 
spectrum obtained from the dipole acceleration with that 
computed only counting the central trajectory of the en- 
semble. Nonetheless, one could now ask the question of 
whether a quantitative agreement between the outcome 
of the TDSE and the results obtained from the Bohmian- 
trajectory computations can be reached. This issue will 
be addressed in the next section. 



C. Ensemble analysis 

In this Section we focus on what distributions of 
Bohmian trajectories lead a quantitative agreement with 
the TDSE solution and the Bohmian-trajectory compu- 
tation. For such a purpose, we have considered a set 
of initial conditions obtained from uniformly generated 
random numbers within a certain interval [— Xc, Xc], sym- 
metric with regard to the origin, a; = 0. This interval has 
been chosen such that the probability density related to 
the initial wavefunction, ^'(a::, 0), is recovered up to over 
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FIG. 6: (Color online) Left panels: Dipole acceleration [panel (a)] and its power spectra [panel (d)] computed from the TDSE 
for a model consisting of a soft-core atomic potential and a sin '^-shaped laser-field pulse with frequency ujq = 0.057 a.u., and 
intensities Eo = 0.05 a.u. (dashed red line) and -Eo = 0.075 a.u. (black solid line). Panels (b) and (c): Bohmian trajectories 
obtained for the same field and atomic parameters as in the left panels, for peak intensities Eo = 0.05 a.u. and Eo = 0.075 
a.u., respectively. Panels (e) and (f): Power spectra from some selected trajectories from among the sets in panels (b) and (c). 
Panel (e) shows the spectra from the trajectories starting at a;(0) = a.u. (black solid line) and x(0) = —3 a.u. (red dashed 
line), for Eo = 0.05 a.u., while panel (f) exhibits spectra computed from the trajectories starting at a;(0) = a.u. (black solid 
line), a;(0) = 1.8 a.u. (blue dotted line) and a;(0) = —3 a.u. (red dashed line), for Eo ~ 0.075 a.u. Colors/line-styles are in 
correspondence with the trajectories displayed at the respective top panels. 



99.8%. For the parameters employed in this work, we 
have verified that all the relevant probability density is 
incorporated if we take Xc — 4.102 a.u. These random 
numbers are multipUed by weights, which are chosen in 
such a way that the probabiUty density related to the 
ground state of the soft-core potential, |V'o(2^)P, is mim- 
icked. For our parameters, this amounts to fourteen in- 
tervals at each side of the origin, x = 0, and around 
N w 1,000 random numbers in total. Each of these 
numbers thus corresponds to the initial condition for a 
Bohmian trajectory. 

This statistical prescription has been used to compute 
the dipole acceleration according to Eq. (|12p for the soft- 
core potential and the trapezoidal pulse. The same field 
strengths as in Fig. [1] have been used. In Figs. El^a) 
and (b), we compare aB(t) with the dipole acceleration 
computed from the TDSE. The outcomes of both ap- 
proaches exhibit an excellent agreement, both qualitative 
and quantitative. This agreement is also observed for the 
corresponding high-order harmonic spectra, displayed in 



panels (c) and (d). 

After having established the good agreement between 
both computations, a legitimate question is what regions 
of the core are the most important for HHG. To investi- 
gate this issue, we have gradually varied the range of the 
aforementioned random-number distribution by changing 
the boundary Xc- These results are displayed in Fig. [5] 
for the time-dependent dipole acceleration and using the 
same intensities as above. As seen, overall, the time 
dependence of all the distributions considered is very 
similar. What is affected is the amplitude of the time- 
dependent acceleration. In fact, if only the acceleration 
along the central trajectory is taken (red lines in the fig- 
ure) , we see that this amplitude is up to five times larger 
than that determined by the acceleration computed with 
the TDSE. The main effect of increasing the range of the 
integration around this central trajectory is to decrease 
this amplitude. Indeed, for both intensities a quantita- 
tive agreement is found for Xc = 2.344 a.u., which cor- 
responds to about 96% of the total probability density. 
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FIG. 7: (Color online) Dipole acceleration for a model con- 
sisting of a soft-core atomic potential Vat{x) and a flat-top 
trapezoidal pulse with frequency ljq = 0.057 a.u and field 
strengths Eo = 0.05 a.u. [panel (a)] and Eq = 0.075 a.u. 
[panel (b)], computed from the TDSE (black solid line) and 
from a Bohmian trajectory simulation employing 987 trajec- 
tories (red dashed line). The spectra corresponding to (a) and 
(b) are represented respectively in panels (c) and (d). 



However, if one is interested in a qualitative comparison 
it is fully legitimate to assume, to first approximation, 
that the core is a source term located around x{0) — 0. 
In fact, this is the key assumption of the existing ana- 
lytic strong-field approaches (e.g., the SFA), which is in 
agreement with our previous results in [30| . 

Finally, in Fig. [S] we investigate how this infiuences 
the high-harmonic spectra. Due to the above-mentioned 
difference in amplitude, the harmonics in the spectrum 
from the central trajectory are several orders of magni- 
tude higher than those from the TDSE [see Figs.[9ja) and 
(b)]. As the spatial range is increased, we observe that a 
reasonable quantitative agreement is reached for the har- 
monics in the cutoff region, already if 75% of the overall 
probability density is considered {xc = 1.172 a.u.). For 
the below-threshold and low-plateau harmonics, a good 
agreement is reached only if 90% of the total probability 
density is taken {xc — 1.8 a.u.). This is still, however, 
less than the requirements in the time-dependent case. 



D. Effective dynamics induced by numerical filters 

Finally, we will address the type of effective dynam- 
ics induced by filters or window functions, which are 
often used to stress some particular aspects from the 
corresponding spectra. For example, a widespread way 
of eliminating the background introduced by the dipole 
length in the harmonic spectra is to apply a Hann (Han- 
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FIG. 8: (Color online) Dipole acceleration obtained from a 
series of sets of Bohmian trajectories randomly distributed 
within the interval [—Xc,Xc] for two values of the laser field: 
(a) Eo = 0.05 a.u. and (b) Eo = 0.075 a.u. Enlargements 
of the behavior of this quantity within one cycle of the pulse 
are given respectively in panels (c) and (d). The different 
values of Xc used in our simulations are labeled with different 
types of colors/line-styles. For details, see the legend in the 
lower panels (in parenthesis, we provide the total number of 
Bohmian trajectories used in each case). In all cases, the 
black sold line represents the result obtained from the TDSE 
and the red dashed line the result from the central Bohmian 
trajectory. 



ning) filter, 



w{t) = 0.5 
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(22) 



This filter is specially useful to force (artificially) the 
time-dependent function g{t) and its derivative to start 
at zero and decay to zero at the end of the simulation in 
a smooth fashion. With it, Eq. ((S]) reads as 



g{t)w{t) e 



(23) 



A legitimate question related to this kind of spectra is 
therefore if they can be associated with a certain dy- 
namical behavior of the system. If g{t) refers to the 
expectation values x{t) or d{t) of the dipole length or 
acceleration, then the effective values Xcff{t) = w{t)x{t) 
and dcff{t) = w{t)d{t) say very little about the behavior 
of an also effective wave function that suits such val- 
ues. In contrast, if g{t) refers to the Bohmian trajecto- 
ries, one can reconstruct an effective dynamics, since (i) 
Xcs{t) — wlt)x{t) is a trajectory (actually, just the prod- 
uct of the actual Bohmian trajectory and the window 
function) and, as seen above, (ii) an ensemble of such 
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FIG. 9: (Color online) Power spectra obtained from the TDSE 
(black solid line) and a series of Bohmian trajectories (red 
dashed line) in which the range of initial conditions have been 
restricted to the intervals [— a;c,a;c], with: Panels (a) and (b): 
Xc = (central trajectory); panels (c) and (d): Xc ~ 1.2 a.u.; 
panels (e) and (f): Xc = 1.8 a.u. In the upper and lower rows 
we show the results for driving-field intensities Eo — 0.05 a.u. 
and Eo = 0.075 a.u., respectively. 



trajectories would allow us to understand how the effec- 
tive wave function (or probability density) would flow in 
the configuration space. That is, it would be like having 
the wave function evolving in an effective potential and 
being affected by an effective pulse. 

In order to show the reconstruction procedure, in 
Fig. [TU] we have considered two trajectories for each of 
the cases shown in Fig. [2] For lower intensities, the tra- 
jectories start and end at their initial position and there- 
fore the window function just stresses this fact (notice 
that wit) produces a reset of the trajectory to zero and 
therefore the effective trajectory has to be shifted so that 
the initial position coincides with that of the respective 
Bohmian trajectory). The effect of w{t) on the Bohmian 
trajectory is just to make it to evolve around the shape 
of this function [see the trajectory starting at x{0) = —3 
a.u., for example, in panel (a)]. In the case of low inten- 
sities and trajectories starting at a::(0) = 0, the spectra 
look essentially the same [see panel (c)]. If the initial po- 
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FIG. 10: (Color online) Bohmian trajectories obtained with 
a flat-top pulse with intensity: (a) Eo = 0.05 a.u. and (b) 
Eo = 0.075 a.u. In panel (a) the trajectories start at x{0) — 
a.u. (black solid line) and a;(0) = —3 a.u. (green dotted 
line); the respective reconstructed effective trajectories are 
denoted with the red dashed line and the blue dash-dotted 
line. In panel (b) the trajectories start at a;(0) = a.u. 
(black solid line) and x(0) = 1.8 a.u. (green dotted line); 
the respective reconstructed effective trajectories are denoted 
with the red dashed line and the blue dash-dotted line. The 
power spectra for each trajectory and its effective homolo- 
gous are given in the four panels below. For the low field 
intensity (£"0 ~ 0.05 a.u.): (c) for x(0) = a.u. and (e) for 
2;(0) — —3 a.u. For the high field intensity {Eo = 0.075 a.u.): 
(d) for x(0) = a.u. and (f) for a:;(0) = 1.8 a.u. Colors/line- 
styles are in correspondence with the trajectories displayed at 
the respective top panels. 



sition is nonvanishing [such as x(0) — —3 a.u., displayed 
in panel (e)], the filter leads to a visible reduction in 
the background. Physically, this happens because of the 
constraint upon the final position Xi{Tf)oi the Bohmian 
trajectory Xi(t) in question imposed by the filter. This 
reduction can be easily spotted in the frequency region 
beyond the cutoff. 

For higher driving-field intensities, however, the fil- 
ter leads to more dramatic effects. In the case of con- 
fined trajectories, the same comments as for low intensi- 
ties hold, as seen in panels (b), regarding the trajectory 
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dynamics, or (d), regarding their spectra (although the 
background is remarkably lower). For the trajectories 
associated with irreversible ionization, though, dynam- 
ics are very different: the Bohmian trajectories goes out 
of the core region, but the effective trajectory obtained if 
the filter is present comes back to its initial position after 
undergoing longer excursions far away. Accordingly, the 
corresponding spectra not only displays more interfer- 
ence features, but the continuous background decreases 
in several orders of magnitude, as seen in panel (f). 



IV. CONCLUSIONS 

Here we have applied Bohmian mechanics to the mod- 
eling of high-order harmonic generation and to the under- 
standing of its local dynamics. Bohmian trajectories have 
been used to monitor the evolution of the wavefunction 
and following that, in the computation of harmonic spec- 
tra. Bohmian trajectories are directly extracted from the 
time-dependent Schrodinger equation. Hence, no physi- 
cal approximations have been made, and that the influ- 
ence of the external field and of the binding potential 
is fully included. Furthermore, they allow a detailed in- 
vestigation of the space-time flow of the time-dependent 
wavefunction, and to associate specific spatial regions of 
this wavefunction to specific harmonic spectra. 

In this work, we have monitored the spectra obtained 
from Bohmian trajectories starting at different spatial re- 
gions, and have that shown that the innermost Bohmian 
trajectory, that is, that starting at x{0) — 0, most closely 
replicates the harmonic spectra on a qualitative level. 
The remaining trajectories mainly contribute to blur up 
these contributions, leading to deviations in the harmonic 
spectra. This issue has been briefiy discussed in our re- 
cent publication [30!|. Here, however, we provide more 
detail and show that this observation holds regardless 
of the type of binding potential, or of the driving laser 
field. This is exemplified by computations for long- and 
short-range potentials in flat-top and sinusoidal pulses. 
Furthermore, in order to provide evidence for the valid- 
ity of our approach, we have shown that a quantitative 
reproduction of the harmonic spectrum obtained from 
the TDSE may be achieved by taking an ensemble of 
Bohmian trajectories whose initial distribution matches 
that of the initial probability density employed in the 
TDSE computation. 

Moreover, one may gain additional insight on how the 
range and the internal structure (e.g., the number of 
bound states and their eigenenergies) of a given potential 
influence ionization dynamics and also HHG by mapping 
these dynamics from the motion of the Bohmian trajec- 
tories. By comparing a long range, one-dimensional soft- 
core potential with truncated models, we observe that 
increased pathways to the continuum, in form of high ly- 
ing bound states, and a reduced potential tail leads to 
more Bohmian trajectories leaving the core region. This 
irreversible ionization causes a strong background on the 



harmonic spectra of the corresponding Bohmian trajec- 
tories, which masks the remaining features. This is in 
agreement with the early work by Krause and co-workers 
|2l| . in which, in the context of the TDSE, such a back- 
ground is attributed to electron density near the spatial 
boundaries inside which the wavefunction is being prop- 
agated. 

This observation has motivated us to revisit the rea- 
sons for the acceleration form of the dipole operator being 
preferable when performing strong field computations in 
this context, especially for stronger driving field intensi- 
ties. We have found that the dipole acceleration effec- 
tively works as a filter that removes from the statistics 
Bohmian trajectories away from the core. These trajec- 
tories lead to the above-mentioned background and are 
over-emphasized in the expectation value of the dipole 
length. Employing Bohmian trajectories, we have shown 
that this effect is similar to that caused by a Banning 
filter, which is commonly used when computing HHG 
spectra from the dipole length. This has been demon- 
strated by back-transforming the spectra of individual 
Bohmian trajectories in which this filter has been em- 
ployed. The Banning filter forces the trajectories related 
to irreversible ionization to return to the core at the end 
of the pulse. After reconstructing the flow of the wave- 
function from an ensemble of Bohmian trajectories, we 
have found that, effectively, a Banning fllter modifies this 
flow and suppresses irreversible ionization. For that rea- 
son, it reduces the background one would obtain if using 
the length form of the dipole operator. 

Finally, it is worth mentioning that one of the main 
results discussed in this paper, namely the fact that the 
trajectory corresponding to the innermost part of the 
time-dependent wavefunction leads to HHG spectra with 
a well-defined plateau and cutoff, seems rather counterin- 
tuitive. This is particularly true in view of the description 
of HHG supported by the three-step model. According 
to this model, HHG as the result of the recombination of 
an electron with its parent ion. Hence, a legitimate ques- 
tion is why a Bohmian trajectory that, spatially, never 
leaves the core, leads to spectra in agreement with such 
predictions. 

This seemingly counterintuitive observation may be 
understood if one takes into account that a Bohmian tra- 
jectory and a classical trajectory are very distinct enti- 
ties. Whereas classically one may identify a bound or 
unbound trajectory by looking at the spatial regions it 
occupies, quantum mechanically one can only talk about 
"bound" and "continuum" in terms of the energy, and 
only if one is considering a time independent problem 
(i.e., eigenstates). If the energy of a classical particle 
is negative (with respect to some relative energy maxi- 
mum), its motion will be bound; if its energy is positive, 
it will be unbound. In contrast, if a quantum system is 
represented by a wave packet, part of it can be inside an 
energy well, part of it can be outside, but both consti- 
tute the same (and unique) wave function. Hence, any 
local alteration in a wavefunction will be felt by the whole 
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wavefunction, i.e., non- locally. Hence, quantum mechan- 
ically, one may have "bound" and "continuum" dynamics 
around x = 0, while classically this is not possible. This 
issue has been addressed in our previous publication [SO^ 
using time-frequency analysis. 
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